{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 栅格数据->矢量数据\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from osgeo import gdal\n",
    "import matplotlib.pyplot as plt\n",
    "from osgeo import ogr,osr\n",
    "import numpy as np\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "ename": "SyntaxError",
     "evalue": "invalid syntax (3541198041.py, line 1)",
     "output_type": "error",
     "traceback": [
      "\u001b[1;36m  Cell \u001b[1;32mIn[2], line 1\u001b[1;36m\u001b[0m\n\u001b[1;33m    path_mndwi=\u001b[0m\n\u001b[1;37m               ^\u001b[0m\n\u001b[1;31mSyntaxError\u001b[0m\u001b[1;31m:\u001b[0m invalid syntax\n"
     ]
    }
   ],
   "source": [
    "path_mndwi=\n",
    "dset=gdal.Open(path_mndwi)\n",
    "mndwi_band=dset.GetRasterBand(1)\n",
    "mndwi_array=mndwi_band.ReadAsArray()\n",
    "print(mndwi_array.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(mndwi_array)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "polygon_path=\n",
    "dset = gdal.Open(path_mndwi)       # 读取路径中的栅格数据\n",
    "mndwi = dset.GetRasterBand(1)    # 这个波段就是最后想要转为矢量的波段，如果是单波段数据的话那就都是1\n",
    "prj = osr.SpatialReference()\n",
    "prj.ImportFromWkt(dset.GetProjection())  # 读取栅格数据的投影信息，用来为后面生成的矢量做准备\n",
    "drv = ogr.GetDriverByName(\"ESRI Shapefile\")\n",
    "polygon = drv.CreateDataSource(polygon_path)  # 创建一个目标文件\n",
    "poly_layer = polygon.CreateLayer('mndwi_polygon', srs=prj, geom_type=ogr.wkbMultiPolygon)  # 对shp文件创建一个图层，定义为多个面类\n",
    "newField = ogr.FieldDefn('DN', ogr.OFTInteger)  # 给目标shp文件添加一个字段，用来存储原始栅格的pixel value\n",
    "poly_layer.CreateField(newField)\n",
    "gdal.Polygonize(srcBand=mndwi, maskBand=None, outLayer=poly_layer, iPixValField=0)   # 核心函数，执行的就是栅格转矢量操作\n",
    "### 删除特定DN值类别要素\n",
    "DN_vec = [1]\n",
    "for fea in poly_layer:\n",
    "    DN_value = fea.GetField('DN')\n",
    "    if DN_value not in DN_vec:\n",
    "        poly_layer.DeleteFeature(fea.GetFID())\n",
    "\n",
    "polygon.SyncToDisk()\n",
    "Polygon = None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 矢量数据->栅格数据\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'path_vec' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32md:\\Lenovo\\Desktop\\云南大学\\开源地理信息系统\\Code\\Chapter5\\20231023.ipynb 单元格 7\u001b[0m line \u001b[0;36m7\n\u001b[0;32m      <a href='vscode-notebook-cell:/d%3A/Lenovo/Desktop/%E4%BA%91%E5%8D%97%E5%A4%A7%E5%AD%A6/%E5%BC%80%E6%BA%90%E5%9C%B0%E7%90%86%E4%BF%A1%E6%81%AF%E7%B3%BB%E7%BB%9F/Code/Chapter5/20231023.ipynb#X10sZmlsZQ%3D%3D?line=3'>4</a>\u001b[0m pixel_size \u001b[39m=\u001b[39m \u001b[39m0.00018\u001b[39m\n\u001b[0;32m      <a href='vscode-notebook-cell:/d%3A/Lenovo/Desktop/%E4%BA%91%E5%8D%97%E5%A4%A7%E5%AD%A6/%E5%BC%80%E6%BA%90%E5%9C%B0%E7%90%86%E4%BF%A1%E6%81%AF%E7%B3%BB%E7%BB%9F/Code/Chapter5/20231023.ipynb#X10sZmlsZQ%3D%3D?line=5'>6</a>\u001b[0m \u001b[39m# Open the data source and read in the extent\u001b[39;00m\n\u001b[1;32m----> <a href='vscode-notebook-cell:/d%3A/Lenovo/Desktop/%E4%BA%91%E5%8D%97%E5%A4%A7%E5%AD%A6/%E5%BC%80%E6%BA%90%E5%9C%B0%E7%90%86%E4%BF%A1%E6%81%AF%E7%B3%BB%E7%BB%9F/Code/Chapter5/20231023.ipynb#X10sZmlsZQ%3D%3D?line=6'>7</a>\u001b[0m dset_vec \u001b[39m=\u001b[39m ogr\u001b[39m.\u001b[39mOpen(path_vec)\n\u001b[0;32m      <a href='vscode-notebook-cell:/d%3A/Lenovo/Desktop/%E4%BA%91%E5%8D%97%E5%A4%A7%E5%AD%A6/%E5%BC%80%E6%BA%90%E5%9C%B0%E7%90%86%E4%BF%A1%E6%81%AF%E7%B3%BB%E7%BB%9F/Code/Chapter5/20231023.ipynb#X10sZmlsZQ%3D%3D?line=7'>8</a>\u001b[0m source_layer \u001b[39m=\u001b[39m dset_vec\u001b[39m.\u001b[39mGetLayer()\n\u001b[0;32m      <a href='vscode-notebook-cell:/d%3A/Lenovo/Desktop/%E4%BA%91%E5%8D%97%E5%A4%A7%E5%AD%A6/%E5%BC%80%E6%BA%90%E5%9C%B0%E7%90%86%E4%BF%A1%E6%81%AF%E7%B3%BB%E7%BB%9F/Code/Chapter5/20231023.ipynb#X10sZmlsZQ%3D%3D?line=8'>9</a>\u001b[0m x_min, x_max, y_min, y_max \u001b[39m=\u001b[39m source_layer\u001b[39m.\u001b[39mGetExtent()\n",
      "\u001b[1;31mNameError\u001b[0m: name 'path_vec' is not defined"
     ]
    }
   ],
   "source": [
    "from osgeo import gdal, ogr\n",
    "\n",
    "# Define pixel_size and NoData value of new raster\n",
    "pixel_size = 0.00018\n",
    "\n",
    "# Open the data source and read in the extent\n",
    "dset_vec = ogr.Open(path_vec)\n",
    "source_layer = dset_vec.GetLayer()\n",
    "proj_vec=source_layer.GetSpatialRef()\n",
    "proj_vec_wkt=proj_vec.ExportToWkt()\n",
    "dset_vec=None\n",
    "x_min, x_max, y_min, y_max = source_layer.GetExtent()\n",
    "\n",
    "geo_trans = (x_min, pixel_size, 0, y_max, 0, -pixel_size)\n",
    "# 定义投影\n",
    "wgs84 = osr.SpatialReference(); wgs84.ImportFromEPSG(4326)\n",
    "wgs84_wkt = wgs84.ExportToWkt()\n",
    "0\n",
    "# Create the destination data source\n",
    "x_size = int((x_max - x_min) / pixel_size)\n",
    "y_size = int((y_max - y_min) / pixel_size)\n",
    "dset_raster = gdal.GetDriverByName('GTiff').Create(path_raster, x_size, y_size, 1, gdal.GDT_Byte)\n",
    "dset_raster.SetGeoTransform(geo_trans)\n",
    "dset_raster.SetProjection(wgs84_wkt)\n",
    "band = dset_raster.GetRasterBand(1)\n",
    "band.SetNoDataValue(-999)\n",
    "# Rasterize\n",
    "gdal.RasterizeLayer(dataset=dset_raster, bands=[1], layer=source_layer, options=[\"ATTRIBUTE=DN\"])\n",
    "# gdal.RasterizeLayer(dataset=dset_raster, bands=[1], layer=source_layer, burn_values=[10])\n",
    "dset_raster = None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.18"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
